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Abstract. We present a solution for the classical univariate rational 
interpolation problem by means of (univariate) subresultants. In the 
case of Cauchy interpolation (interpolation without multiplicities), we 
recover explicit formulas for the solution in terms of symmetric functions 
of the input data, hence generalizing the well-known formulas for La- 
grange interpolation. In the case of the osculatory rational interpolation 
(interpolation with multiplicities), we get determinantal expressions in 
terms of the input data. 



1. Cauchy and osculatory interpolation 

The Cauchy interpolation problem or rational interpolation problem con- 
sidered already in [Caul84l| IRos 18451 IPredl9 53] is the following: 

Let K be a field, and o, b G No- Set I = a + b. Given a 
set {xo, . . . , xi\ of I + 1 points in K, and yo, . . . ,y# G K, 
determine -if possible- polynomials A, B G K[x] such that 
deg(^4) < a, deg(B) < b and 

(1) -g (a*) = yi, < i < £. 

This might be considered as a generalization of the classical Lagrange 
interpolation problem for polynomials, where 6 = and a = £. In contrast 
with that case, there is not always a solution to this problem, since for 
instance by setting yo = ■ ■ ■ = y a = 0, the numerator A is forced to be 
identically zero, and therefore the remaining y a +k, 1 < k < £ — a, have to 
be zero as well. 

The obvious generalization of the Cauchy interpolation problem receives 
the name osculatory rational interpolation problem or rational Hermite in- 
terpolation problem: 

Let K be a field, and a, b G No- Set £ = a + b. Given a set 
{xo, . . . , Xk} of k + 1 points in K , oq, . . . , G N such that 
«0 + ' ' ■+cik = ^+1, and j/jj G K, < i < k , < j < ai, not 
all of them equal to zero, determine -if possible- polynomials 
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A, B € K[x] such that deg(A) < a, deg(B) < b and 
/ A\ (J) 

(2) f - J (xi) = jlyy, < i < fc, < j < at. 

This problem has also been extensively studied from both an algorithmic 
and theoretical point of view, see for instance [Sall962, Kahl969, Wuyl975 



TF2000J and the references therein. In all these previous works, the solution 
of the problem comes in some kind of recursive continuous fraction, and 
it is not actually explicit. A unified framework, called rational function 
reconstruction, is presented in the book [vzGG2003] for both the Cauchy 
interpolation and the oscillatory rational interpolation problems. It uses 
remainder sequences in the Extended Euclidean Algorithm. 

Yet an explicit general formula in terms of the input data is not found 
in the literature. In the classical (Cauchy) case, an explicit formula can be 
derived from the results on symmetric operators in a suitable ring of poly- 
nomials presented in [Las2003j, as shown in [Las]. In this paper, we get 
the explicit formulas for the Cauchy interpolants by expressing the subre- 
sultants in terms of the Sylvester sums introduced by Sylvester in |Syll853|. 
Furthermore, we give determinantal expressions for the solution of the os- 
culatory rational interpolation problem in terms of the input data. These 
latter expressions are given as quotients of determinants of matrices. Note 
that there are also explicit formulas for the Hermite interpolants in terms of 
the input data, see for instance [DKS2012J for a presentation of them. But 
yet we do not have a closed formula for subresultants in roots in the multiple 
case, so a generalization of Theorem 11.11 is yet not a forward application of 
Sylvester's multiple sums, and some more work on the subject must be done 
in order to shed light to the general problem. 

We first present the result for the Cauchy interpolation problem. For 
U,VCK, we set R(U, V) := Uueu,vev (« " «)• 

Theorem 1.1. Given (a, b), X := {xo, . . . , xg} and yo, . . . , yi as above. Let 
d be maximal such that < a < d and 

A:= Yl R(^ X ')( II yj)/R(X\X',X') EK[x] 

X'cX,\X'\=d 

is not identically zero. Then a solution for the Cauchy interpolation problem 
(HJ exists if and only if for 

B := ]T R(X",x)( J] yj )/R{X",X\X") £ K[x], 
x"cx,\x"\=e-d xjex" 

one has B(x{) ^ for all < i < t. In that case the (essentially unique) 
solution is given by A/B. 

In particular, when a = £, Theorem 11.11 will specialize to the well-known 

Lagrange interpolation polynomial, associated to the data {(xi, yi)}o<i<£, 
which is equal to 

Th#fr ~ x i) = sr R ( x i x \i x i}) 

kJ%^-^ i- <e yi R( Xi ,x\{ Xi }y 
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Before stating our main result for the osculatory rational interpolation 
problem, we need to set a notation: 

Notation 1.2. Set a, b 6 N such that a + b = I, oq, ■ • • , cife € N such that 
ao + ' ' ' + &k = i + 1; as in ([2]). M^e define 

• X := ((xq, oo); . . . ; afc)) an array of pairs in K X N and Y := 
(yijj < i < /c, < j < a,j\. We call these the input data for the 
osculatory rational interpolation problem. 

• Set u £ N. The generalized Vandermonde or confluent matrix (c.f. 
|Kall9 84] ) of size u + 1 associated to X is the (non-necessarily 
square) matrix V u +i(X) € f^(u+i)x(l+l) defined by 



K+iPO 




u+l 



where for any t, V u +i(xi,t + 1) € K( u + 1 ) x ( t + 1 ) is defined by 

t+i 





1 















1 








V u+ i(xi,t+l) := 


*? 




1 









uxr 1 




• or* 



u+l 



We define the matrix U U+1 {X,Y) £ K(«+ 1 )x(^+ 1 ) associated to X 
and Y as: 



U U+1 (X,Y) := 




E4+i(£fc;2/fc,o,---,2/fc,a fc -i) 



u+l 



where for any t, U u+ \{xi,yifl, . . . ,y% t t) £ if ( u+1 ) x (* +1 ) is defined by 



t+l 


Ui,0 






XiVifi 


XiVi,i + Vi,o 


x i y l . t + Vi,t-i 




xfVi,! +uxf~ 2 y lfi ... 


2jj'=0 (j) X i 3 Vi,t—0 



u+l 



with the convention that when u < j, 3 = 0. 

The next result expresses the solution of the osculatory rational interpo- 
lation problem in terms of the input data as follows: 
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Theorem 1.3. Under the notation above, let d be maximal such that < 
d < a and 



e+i 


1 






1 




V d+1 (X) 




A:= — det 




x d 




U e _ d+1 (X,Y) 






d+l 



t-d+l 



€ K[x] 



is not identically zero. Then a solution for the osculatory rational interpo- 
lation problem ([2]) exists if and only if for 



e+i 



B := det 



V d+l (X) 





U e . d+1 (X,Y) 


1 

x e-d 



d+l 



t-d+l 



€ K[x] 



one has B{xi) ^ for all < i < k. In that case the (essentially unique) 
solution is given by A/B. 

2. Proof of the results 

Let us start by showing that a solution for the rational interpolation 
problem is essentially unique. 

Proposition 2.1. If the osculatory rational interpolation problem ([2|) has 
a solution, then there is a unique solution A/B with gcd(A, B) = 1 and A 
monic. 

Proof. If there is a solution, then, cleaning common factors and dividing by 
the leading coefficient of A, there is a solution satisfying the same degree 
bounds with gcd(A, B) = 1 and A monic. Assume A\jB\ and A2/B2 are 
both solutions of the same type. Then, [A\/ Bx)^'{xi) = {A2/ B2)^\xi) 
implies 

(^i^s - A 2 B 1 /B 1 B 2 ) ij) (x l ) = for < i < k, < j < a h 

which inductively implies that {A\B 2 — A 2 Bi)^ (xi) = for the £ + 1 condi- 
tions. But A\B 2 — A 2 B\ is a polynomial of degree at most £, and therefore 
A\B 2 = A 2 B\. Therefore, A\ = cA 2 and B\ = cB 2 with c € K \ {0}, as 
both A\ and A 2 are monic, then c = 1 and the claim follows. □ 

To prove our results, we first connect the solutions of the osculatory ra- 
tional interpolation problem (and the Cauchy interpolation problem as a 
special case) to subresultants of the following two polynomials: 

• / := ni=o( x ~ x j) aj > which we write / = ^=0 fax 1 , where f i+1 = 1. 
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• g € K[x] the unique polynomial of degree less than or equal to £ 
such that 

9 b) i x i) = 3- Vi,j, < i < k, 0<j<a i} 

which we expand as g = J2i=o 9i x ' \ where gi = for deg(g) < i < £. 
This polynomial is the so-called Hermite interpolation polynomial 
associated to the data 

{(xq, a ), • • • , (x k , a k )}, {yi,j}o<i<k,o<j< ai - 

By assumption j / since we assumed that some of the y^j is 
non-zero. 

For d < £, consider the d-th subresultant polynomial Sresd(f,g) of / and g, 
defined as 

2l+l-2d 



t-d 



e+i-d 



Note that the previous definition makes sense even if deg(g) = m < £, 
and agrees for d < m with the usual definition of subresultant of / and g 
given by the matrix of the right size £ + 1 + m — 2d, since / is monic. For 
m < d < £ we have, according to the definition above, that Sresrf(/, g) = 0, 
and for d = £, Sies e (f,g) =g = Sres m (/,#). 

We have the universal subresultant Bezout identity 

(3) Sves d (f,g)=F d f + G d g, 

where 



t-d 



e+i-d 



and 

2£+l-2d 



t-d 



t+l-d 



Sres d (f,g) := dct 



fi+l ■ 




fd+l-(t-d-l) 






fe+i ■ 


fd+1 


x°f(x) 






9d+i-(e-d) 


x t ~ a g{x) 




9i ■ 


9d+l 


x°g(x) 



F d det 



ft+l ■ 




Jd+l-(t-d-l) 


x t-d-l 




fe+i ■ 


fd+1 


x° 


91 ■ 




9d+l-{l-d) 







9t ■ 


9d+l 






G d := dct 



fi+l ■ 




fd+l-{t-d-l) 







ft+i ■ 


fd+1 





9i 




9d+l-(t-d) 


x e- d 




9i ■ 


9d+l 


x° 
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Observe that deg(G d ) < t — d when G d ^ 0. 

The result below, strongly related to [vzGG2003, Theorem 5.16], expresses 
the existence and uniqueness of the solution of the osculatory rational inter- 
polation problem in terms of the subresultant sequence of / and g. 

Proposition 2.2. Let (a, b) be given, and set d < a to be maximal such 
that Sres d (f,g) ^ 0. Then the osculatory rational interpolation problem ([2]) 
has a solution if and only if Gd(xi) ^ for all < i < k. In this case, the 
(essentially unique) solution is given by 

A = Sres d (f,g) 
B G d ' 

and moreover gcd (Sres^/, g), Gd) = 1. 

Proof. We use [vzGG20031 Th. 5.16]: Problem ([2]) has a solution if and only 
if the minimal row rj = Sjf + tjg in the Extended Euclidean Algorithm such 
that dj := deg(rj) < a satisfies that gcd(rj,tj) = 1, in which case rj/tj is 
the canonical solution. Note that for di := deg(rj), we have dj < a < cL-i 
by the definition of r,- . 

Then, by the Fundamental Theorem of PRS (c.f. [ Colll9671 IBT1971j or 
[GCL1996| Th.7.4]), Sres^ (/, g) and S?es dj _ 1 -i{f , g) are (non-zero) constant 
multiples of rj, and Sres^ (/, g) = for dj < d! < — 1. In any case, if 
d < a is the largest such that Res^/, g) ^ 0, one has Resrf(/, g) = crj for 
some c ^ in K. 

But Res d(/,g) = F d f + G d g with deg(Res d (/, g)) + deg(G d ) < deg(/) im- 
plies by [vzGG2003, Lemma 5.15] that F d = csj and G d = ctj too. More- 
over, by |vzGG2003~| Lemma 3.15 (v)], gcd(F d ,G d ) = 1, which implies that 
gcd (Sies d ( f,g),G d ) = 1 & G d [xi) ^ for all < i < k, by the def- 
inition of /. Therefore, there is a solution to the problem if and only if 
G d (xi) 7^ 0, for all < i < k, and the essentially unique solution is given by 
Sres d (f,g)/G d . ' □ 



Proposition 12.21 has the advantage that it can be applied to produce ex- 
plicit formulae for the rational interpolation problem in terms of the input 
data. We first prove our result on the Cauchy interpolation problem, when 
/ has simple roots. 



Proof of Theorem 11.11 Set X = {xq,...,x{\ for the set of roots of 

/ = ni=o( x — x «) ana - ^ f° r * ne se ^ °f m ro °t s of 5S t ne unique polyno- 
mial of degree bounded by £ which satisfies g(xi) = y. L for < i < £. 
Using Proposition 12.21 above, it is sufficient to prove that Sres d (f, g) = A 
and G d = B. We use Sylvester's single-sum formula in roots for Sres^ (see for 
instance the original paper of Sylvester |Syll853 Art. 21] or the many other 
references on the topic), and also for G d ( |Syll853j Art. 29], or |KS2012| . 
Remark after Lemma 6): 
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Sres d (f,g) = ^ R(x,X 



\X'\=d 



„ R(X\X',Z) 
R(X\X',X')' 



V R(x X') 9 ^ Xi) 
L> [ ' ' R(X\X',X>Y 

\X'\=d \ \ i J 



G d = {-lf- d Yl R ^ X " 

\X"\=l-d 



R(X", Z) 



R[X",X\X n ) 
R(X", Z) 



\x"\=e-d 



R(X",X\X") 



- ^ n ^ ,X) R(X",X\X"Y 
\X"\=t-d \ i \ J 

where both X', X" C X. We get that A = Sres d (f,g) and B = G d if we 
substitute = yj for all i above. □ 

Observe that in the previous theorem, if a = £ then Siesi(f,g) = j / 
and the solution is given by the Lagrange interpolation polynomial, as 
mentioned in the introduction. 

To prove Theorem 11.31 for the osculatory rational interpolation problem 
we need the following notation, that we already used in [DKS2012J. 

Notation 2.3. Recall X = ((xo, oq); . . . ; (x^, a&)) . Given a polynomial h(z) 
and u € N, the generalized Wronskian of size u + 1 associated to X is the 
(non-necessarily square) matrix Wh,u+l(X) £ i^( u + 1 ) x (^+ 1 ) defined by 



e+i 



W Ku+1 {X) 



u+1 , 



where for any t, Wh u +\(xi,t + 1) G i^("+ 1 ) x (*+ 1 ) is defined by 



t+i 





h(xi) 


h'(xi) 


t\ 


W h , u +l{xi,t+l) := 


[zh)(xi) 


{zh)'(Xi) . 






(z u h)(Xi) 


(z U h)'(Xi) . 


(z"/*) (t) (*,) 



with the convention that when k < j, Q)x* J = 0. 
First we prove the following lemma: 

Lemma 2.4. Using the notation above, if there is a solution for the oscu- 
latory rational interpolation problem, then the essentially unique solution is 



s 



CARLOS D' ANDREA, TERESA KRICK, AND AGNES SZANTO 



given by 



1+1 



A 
B 







1 




V d+ i{X) 




det 




x d 




W 9tt - d+x (X) 






e-d+i 



e+i 



det 







V d+1 (X) 







1 


W 9tt - d+x (X) 






x e '~ d 



d+l 



-d+1 



where d < a is maximum such that the numerator does not vanish. 

Proof. We apply Proposition l2.2l for the maximum d < a such that Sresd(/, g) ^ 
0. In [DKS2012, Theorem 2.5] we proved that 



Sres d (/,<?) = {-lf+ x - d det (V e+1 (X)) 1 det 



t+l 


1 






1 




V d +i{X) 




det 




x d 




W g ^ d +i(X) 






d+l 



-d+l 



To prove a similar statement about consider the following matrices: 

2t-d+2 



M f := 



2£-d+2 


fo ■ 


■ h+i 









h ... 


ft+i 






t-d . M 9 ■■ = 



go ■ 


■ gi 








go ■ 


■ ge 


x e-d 



and 



2£-d+2 



u d 



I d+l 



M f 



d+l 

l-d 

e+i-d 



where Id+i is the (d+l) x (2£ — d + 2) matrix with the identity matrix on 
the left and zero otherwise. Then from the definition of Gd we have that 



G d = det(U d ). 
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Also, similarly as in [DKS2012,, Theorem 2.5], we have 



e-d 1 



2i-d+2 



d+1 

e-d 

t+l—d 



M f 



t+l 


e-d+i 




V d +i{X) 


* 





V-U-d+l{X) 





t+l 





M> 


















t-d+1 






1 








W gM ^ d (X) 


* 


x l-d 



d+1 



e-d 



e+i-d 



where Mj- is a triangular matrix with fg + \ = 1 in its diagonal. This implies 
that 

e+i i 



G d = (-iy- d dct(V e+1 (X)) Met 



V d +i(X) 





W g ^ d +i(X) 


1 

x t-d 



d+1 



-1-d 



Now we can easily prove Theorem 11.3 



□ 



Proof of Theorem II. 3L We simply compute in the previous lemma the 
entries of the matrix W gi i-d+i(X), applying Leibniz rule and the fact that 
g^ j \xi) = (t - j)\ yij for < i < k and < j < af. 

(z u g)( t \x l ) A (u 



t\ 



U—J 

j=0 V ' 



□ 
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plained us part of the results in [Las2003j. 



References 

[BT1971] Brown, W. S.; Traub, J. F. On Euclid's algorithm and the theory of subresul- 
tants. J. Assoc. Comput. Mach. 18 (1971), 505-514. 

[Caul841] A. L. Cauchy. Memoire sur les fonctions alternees et les sommes alternees. 
Exercices d'analyse et de phys. math. (1841) 151-159. 

[ColU967] George E. Collins. Subresultants and Reduced Polynomial Remainder Se- 
quences. J. ACM 14, 1 (1967), 128-142. 

[DKS2012] D'Andrea, Carlos; Krick, Teresa; Szanto, Agnes. Subresultants in multiple 
roots. Linear Algebra and its Applications, to appear. 

[GCL1996] Geddes, Keith; Czapor, S.; Labahn, G. Algorithms for Computer Algebra. 
Kluwer Academic Publishers, 1996. 

[Kahl969] Kahng, S. W. Osculatory interpolation. Math. Comp. 23 1969 621-629. 



10 



CARLOS D' ANDREA, TERESA KRICK, AND AGNES SZANTO 



[Kall984] Kalman, D. The generalized Vandermonde matrix. Math. Mag. 57 (1984), no. 
1, 15-21. 

[KS2012] Krick, Teresa; Szanto, Agnes. Sylvester's double sums: an inductive proof of 
the general case. J. Symbolic Comput. 47 (2012), no. 8, 942-953. 

[Las2003] Lascoux, Alain. Symmetric functions and combinatorial operators on polyno- 
mials. CBMS Regional Conference Series in Mathematics, 99. Published for 
the Conference Board of the Mathematical Sciences, Washington, DC; by the 
American Mathematical Society, Providence, RI, 2003. 

[Las] Lascoux, Alain. Notes on Interpolation in one and several variables. 

|http : //igm . univ-mlv . fr/~al /ARTICLES/int erp . dvi . gz | 

[Predl953] Predonzan, Arno. Su una formula d'interpolazione per le funzioni razionali. 
Rend. Sem. Mat. Univ. Padova 22, (1953). 417-425. 

[Rosl845] G. Rosenhain. Neue Darstellung der Resultante der Elimination von z aus 
zwei algebraische Gleichungen. Crelle J. 30 (1845) 157-165. 

[Sall962] Salzer, Herbert E. Note on osculatory rational interpolation. Math. Comp. 16 
1962 486-491. 

[Syll853] Sylvester, James Joseph. On a theory of syzygetic relations of two rational 
integral functions, comprising an application to the theory of Sturm 's function 
and that of the greatest algebraical common measure. Philosophical Transac- 
tions of the Royal Society of London, Part III (1853), 407-548. Appears also 
in Collected Mathematical Papers of James Joseph Sylvester, Vol. 1, Chelsea 
Publishing Co. (1973) 429-586. 

[TF2000] Tan, Jieqing; Fang, Yi. Newton-Thtele's rational interpolants. Computational 
methods from rational approximation theory (Wilrijk, 1999). Numer. Algo- 
rithms 24 (2000), no. 1-2, 141-157. 

[vzGG2003] von zur Gathen, Joachim; Gerhard, Jiirgen. Modern computer algebra. Second 
edition. Cambridge University Press, Cambridge, 2003. 

[Wuyl975] Wuytack, Luc. On the osculatory rational interpolation problem. Math. Com- 
put. 29 (1975), 837-843. 

Universitat de Barcelona, Departament d'Algebra i Geometria. Gran Via 
585, 08007 Barcelona, Spain. 
E-mail address: cdandrea@ub.edu 
URL: http: //atlas .mat .ub . es/personals/dandrea 

Departamento de Matematica, Facultad de Ciencias Exactas y Naturales, 
Universidad de Buenos Aires and IMAS, CONICET, Argentina 
E-mail address: krick@dm.uba.ar 
URL: http://mate.dm.uba.ar/~krick 

Department of Mathematics, North Carolina State University, Raleigh, 
NC 27695 USA 

E-mail address: aszanto@ncsu.edu 
URL: www4.ncsu. edu/~ aszanto 



